Introduction

In this document, we will outline the Bayesian analogs of the statistical analyses described here (and in the previous course lecture).

Helper functions.

# this function extracts results from different models and generate results of the same format to be used in visualizations
tidy.wrapper = function(model) {
  if (class(model) == "lm") {
    tidy(model, conf.int = TRUE) %>%
      select(-c(statistic, p.value)) %>%
      mutate(model = "Frequentist") %>%
      select(model, everything())
  } else if (class(model) == "brmsfit") {
    tidy(model) %>%
      filter(effect == "fixed") %>%
      select(c(term:conf.high)) %>%
      mutate(model = "Bayesian") %>%
      select(model, everything())
  } else {
    stop("unknown model class")
  }
}

TODO: explain the dataset here

Dataset 1

Let’s first load and look at the data in a table

# Fumemng -> Abhraneel I'm using the Rproj in the root folder, although my path is identical...

dataset = readr::read_csv("02-bayesian_stats/data/blinded.csv")
## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   experiment = col_double(),
##   condition = col_character(),
##   effectiveness = col_double()
## )
# dataset = readr::read_csv("~/Documents/Github/chi22-course/data/blinded.csv")

kable(head(dataset))
experiment condition effectiveness
1 no_graph 7
1 graph 6
1 no_graph 9
1 no_graph 4
1 graph 6
1 no_graph 7
dataset %>% 
  mutate(effectiveness = fct_rev(factor(effectiveness, levels = 1:9))) %>%
  
  # stacked bar plot
  ggplot(aes(x = condition, fill = effectiveness)) +
  geom_bar(position = "stack", stat="count") +
  
  # plot data for different experiments as small multiples
  facet_wrap( ~ experiment) +
  
  # grey color scale is robust for colorblind
  scale_fill_brewer(palette="Greys", drop = FALSE) +
  
  # horizontal plot
  coord_flip() +
  
  # legend
  guides(fill = guide_legend(reverse = TRUE)) 

As we can see above, the original dataset contains results from four different experiments. For the purposes of this lecture, we will confine ourselves to the first experiment.

exp1.data = dataset %>%
  filter(experiment == 1)

exp1.data
## # A tibble: 123 × 3
##    experiment condition effectiveness
##         <dbl> <chr>             <dbl>
##  1          1 no_graph              7
##  2          1 graph                 6
##  3          1 no_graph              9
##  4          1 no_graph              4
##  5          1 graph                 6
##  6          1 no_graph              7
##  7          1 graph                 7
##  8          1 no_graph              7
##  9          1 no_graph              7
## 10          1 graph                 6
## # … with 113 more rows

Model 1. Wilcoxon signed rank test

This is a non-parametric test which we will skip for now. Although, there exists Bayesian non-parametric methods, they are more advanced for this lecture.

Model 2. t-test

This is the linear model equivalent for the paired sample t-test

dataset1.lm.freqt <-
  lm(
    effectiveness ~ condition - 1, 
    data = exp1.data
  )
  exp1.data %>%
     ggplot(., aes(x = effectiveness)) +
  geom_histogram(fill = '#351c75', color = NA)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

dataset1.brm.bayesiant.formula <- bf(
    effectiveness ~ condition,
    family = student()
  )
kable(as.tibble(get_prior(dataset1.brm.bayesiant.formula, data = exp1.data)))
## Warning: `as.tibble()` was deprecated in tibble 2.0.0.
## Please use `as_tibble()` instead.
## The signature and semantics have changed, see `?as_tibble`.
prior class coef group resp dpar nlpar bound source
b default
b conditionno_graph default
student_t(3, 7, 2.5) Intercept default
gamma(2, 0.1) nu default
student_t(3, 0, 2.5) sigma default
dataset1.brm.bayesiant.priors = c(
      prior(normal(0, 1), class = "b"), # there's a lot of data so even fairly "strong" priors are going to not matter so much here
      prior(student_t(3, 0, 1), class = "sigma")
      )
dataset1.brm.bayesiant.priorchecks <-
  brm(
    dataset1.brm.bayesiant.formula,
    prior = dataset1.brm.bayesiant.priors,
    data = exp1.data,
    backend = "cmdstanr",
    sample_prior = "only",
    file = "02-bayesian_stats/rds/dataset1.brm.bayesiant.priorchecks.rds"
  )
## Compiling Stan program...
## Start sampling
## Running MCMC with 4 sequential chains...
## 
## Chain 1 Iteration:    1 / 2000 [  0%]  (Warmup) 
## Chain 1 Iteration:  100 / 2000 [  5%]  (Warmup) 
## Chain 1 Iteration:  200 / 2000 [ 10%]  (Warmup) 
## Chain 1 Iteration:  300 / 2000 [ 15%]  (Warmup) 
## Chain 1 Iteration:  400 / 2000 [ 20%]  (Warmup) 
## Chain 1 Iteration:  500 / 2000 [ 25%]  (Warmup) 
## Chain 1 Iteration:  600 / 2000 [ 30%]  (Warmup) 
## Chain 1 Iteration:  700 / 2000 [ 35%]  (Warmup) 
## Chain 1 Iteration:  800 / 2000 [ 40%]  (Warmup) 
## Chain 1 Iteration:  900 / 2000 [ 45%]  (Warmup) 
## Chain 1 Iteration: 1000 / 2000 [ 50%]  (Warmup) 
## Chain 1 Iteration: 1001 / 2000 [ 50%]  (Sampling) 
## Chain 1 Iteration: 1100 / 2000 [ 55%]  (Sampling) 
## Chain 1 Iteration: 1200 / 2000 [ 60%]  (Sampling) 
## Chain 1 Iteration: 1300 / 2000 [ 65%]  (Sampling) 
## Chain 1 Iteration: 1400 / 2000 [ 70%]  (Sampling) 
## Chain 1 Iteration: 1500 / 2000 [ 75%]  (Sampling) 
## Chain 1 Iteration: 1600 / 2000 [ 80%]  (Sampling) 
## Chain 1 Iteration: 1700 / 2000 [ 85%]  (Sampling) 
## Chain 1 Iteration: 1800 / 2000 [ 90%]  (Sampling) 
## Chain 1 Iteration: 1900 / 2000 [ 95%]  (Sampling) 
## Chain 1 Iteration: 2000 / 2000 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 Iteration:    1 / 2000 [  0%]  (Warmup) 
## Chain 2 Iteration:  100 / 2000 [  5%]  (Warmup) 
## Chain 2 Iteration:  200 / 2000 [ 10%]  (Warmup) 
## Chain 2 Iteration:  300 / 2000 [ 15%]  (Warmup) 
## Chain 2 Iteration:  400 / 2000 [ 20%]  (Warmup) 
## Chain 2 Iteration:  500 / 2000 [ 25%]  (Warmup) 
## Chain 2 Iteration:  600 / 2000 [ 30%]  (Warmup) 
## Chain 2 Iteration:  700 / 2000 [ 35%]  (Warmup) 
## Chain 2 Iteration:  800 / 2000 [ 40%]  (Warmup) 
## Chain 2 Iteration:  900 / 2000 [ 45%]  (Warmup) 
## Chain 2 Iteration: 1000 / 2000 [ 50%]  (Warmup) 
## Chain 2 Iteration: 1001 / 2000 [ 50%]  (Sampling) 
## Chain 2 Iteration: 1100 / 2000 [ 55%]  (Sampling) 
## Chain 2 Iteration: 1200 / 2000 [ 60%]  (Sampling) 
## Chain 2 Iteration: 1300 / 2000 [ 65%]  (Sampling) 
## Chain 2 Iteration: 1400 / 2000 [ 70%]  (Sampling) 
## Chain 2 Iteration: 1500 / 2000 [ 75%]  (Sampling) 
## Chain 2 Iteration: 1600 / 2000 [ 80%]  (Sampling) 
## Chain 2 Iteration: 1700 / 2000 [ 85%]  (Sampling) 
## Chain 2 Iteration: 1800 / 2000 [ 90%]  (Sampling) 
## Chain 2 Iteration: 1900 / 2000 [ 95%]  (Sampling) 
## Chain 2 Iteration: 2000 / 2000 [100%]  (Sampling) 
## Chain 2 finished in 0.0 seconds.
## Chain 3 Iteration:    1 / 2000 [  0%]  (Warmup) 
## Chain 3 Iteration:  100 / 2000 [  5%]  (Warmup) 
## Chain 3 Iteration:  200 / 2000 [ 10%]  (Warmup) 
## Chain 3 Iteration:  300 / 2000 [ 15%]  (Warmup) 
## Chain 3 Iteration:  400 / 2000 [ 20%]  (Warmup) 
## Chain 3 Iteration:  500 / 2000 [ 25%]  (Warmup) 
## Chain 3 Iteration:  600 / 2000 [ 30%]  (Warmup) 
## Chain 3 Iteration:  700 / 2000 [ 35%]  (Warmup) 
## Chain 3 Iteration:  800 / 2000 [ 40%]  (Warmup) 
## Chain 3 Iteration:  900 / 2000 [ 45%]  (Warmup) 
## Chain 3 Iteration: 1000 / 2000 [ 50%]  (Warmup) 
## Chain 3 Iteration: 1001 / 2000 [ 50%]  (Sampling) 
## Chain 3 Iteration: 1100 / 2000 [ 55%]  (Sampling) 
## Chain 3 Iteration: 1200 / 2000 [ 60%]  (Sampling) 
## Chain 3 Iteration: 1300 / 2000 [ 65%]  (Sampling) 
## Chain 3 Iteration: 1400 / 2000 [ 70%]  (Sampling) 
## Chain 3 Iteration: 1500 / 2000 [ 75%]  (Sampling) 
## Chain 3 Iteration: 1600 / 2000 [ 80%]  (Sampling) 
## Chain 3 Iteration: 1700 / 2000 [ 85%]  (Sampling) 
## Chain 3 Iteration: 1800 / 2000 [ 90%]  (Sampling) 
## Chain 3 Iteration: 1900 / 2000 [ 95%]  (Sampling) 
## Chain 3 Iteration: 2000 / 2000 [100%]  (Sampling) 
## Chain 3 finished in 0.0 seconds.
## Chain 4 Iteration:    1 / 2000 [  0%]  (Warmup) 
## Chain 4 Iteration:  100 / 2000 [  5%]  (Warmup) 
## Chain 4 Iteration:  200 / 2000 [ 10%]  (Warmup) 
## Chain 4 Iteration:  300 / 2000 [ 15%]  (Warmup) 
## Chain 4 Iteration:  400 / 2000 [ 20%]  (Warmup) 
## Chain 4 Iteration:  500 / 2000 [ 25%]  (Warmup) 
## Chain 4 Iteration:  600 / 2000 [ 30%]  (Warmup) 
## Chain 4 Iteration:  700 / 2000 [ 35%]  (Warmup) 
## Chain 4 Iteration:  800 / 2000 [ 40%]  (Warmup) 
## Chain 4 Iteration:  900 / 2000 [ 45%]  (Warmup) 
## Chain 4 Iteration: 1000 / 2000 [ 50%]  (Warmup) 
## Chain 4 Iteration: 1001 / 2000 [ 50%]  (Sampling) 
## Chain 4 Iteration: 1100 / 2000 [ 55%]  (Sampling) 
## Chain 4 Iteration: 1200 / 2000 [ 60%]  (Sampling) 
## Chain 4 Iteration: 1300 / 2000 [ 65%]  (Sampling) 
## Chain 4 Iteration: 1400 / 2000 [ 70%]  (Sampling) 
## Chain 4 Iteration: 1500 / 2000 [ 75%]  (Sampling) 
## Chain 4 Iteration: 1600 / 2000 [ 80%]  (Sampling) 
## Chain 4 Iteration: 1700 / 2000 [ 85%]  (Sampling) 
## Chain 4 Iteration: 1800 / 2000 [ 90%]  (Sampling) 
## Chain 4 Iteration: 1900 / 2000 [ 95%]  (Sampling) 
## Chain 4 Iteration: 2000 / 2000 [100%]  (Sampling) 
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 1.0 seconds.
dataset1.bayesiant.yprior <-
  posterior_predict(dataset1.brm.bayesiant.priorchecks)

ggplot() +
  geom_density(aes(x = dataset1.bayesiant.yprior),
               color = '#351c75',
               alpha = .5,
               size = 1) +
  xlab('prior draws') +
  ggtitle('prior preditive check')

dataset1.brm.bayesiant = 
  brm(dataset1.brm.bayesiant.formula,
    prior = dataset1.brm.bayesiant.priors,
    data = exp1.data,
    #You may not need this if rstan works with your computer
    backend = "cmdstanr", 
    #Save the model
    file = "02-bayesian_stats/rds/dataset1.brm.bayesiant.rds"
  )
## Compiling Stan program...
## Start sampling
## Running MCMC with 4 sequential chains...
## 
## Chain 1 Iteration:    1 / 2000 [  0%]  (Warmup) 
## Chain 1 Iteration:  100 / 2000 [  5%]  (Warmup) 
## Chain 1 Iteration:  200 / 2000 [ 10%]  (Warmup) 
## Chain 1 Iteration:  300 / 2000 [ 15%]  (Warmup) 
## Chain 1 Iteration:  400 / 2000 [ 20%]  (Warmup) 
## Chain 1 Iteration:  500 / 2000 [ 25%]  (Warmup) 
## Chain 1 Iteration:  600 / 2000 [ 30%]  (Warmup) 
## Chain 1 Iteration:  700 / 2000 [ 35%]  (Warmup) 
## Chain 1 Iteration:  800 / 2000 [ 40%]  (Warmup) 
## Chain 1 Iteration:  900 / 2000 [ 45%]  (Warmup) 
## Chain 1 Iteration: 1000 / 2000 [ 50%]  (Warmup) 
## Chain 1 Iteration: 1001 / 2000 [ 50%]  (Sampling) 
## Chain 1 Iteration: 1100 / 2000 [ 55%]  (Sampling) 
## Chain 1 Iteration: 1200 / 2000 [ 60%]  (Sampling) 
## Chain 1 Iteration: 1300 / 2000 [ 65%]  (Sampling) 
## Chain 1 Iteration: 1400 / 2000 [ 70%]  (Sampling) 
## Chain 1 Iteration: 1500 / 2000 [ 75%]  (Sampling) 
## Chain 1 Iteration: 1600 / 2000 [ 80%]  (Sampling) 
## Chain 1 Iteration: 1700 / 2000 [ 85%]  (Sampling) 
## Chain 1 Iteration: 1800 / 2000 [ 90%]  (Sampling) 
## Chain 1 Iteration: 1900 / 2000 [ 95%]  (Sampling) 
## Chain 1 Iteration: 2000 / 2000 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 Iteration:    1 / 2000 [  0%]  (Warmup) 
## Chain 2 Iteration:  100 / 2000 [  5%]  (Warmup) 
## Chain 2 Iteration:  200 / 2000 [ 10%]  (Warmup) 
## Chain 2 Iteration:  300 / 2000 [ 15%]  (Warmup) 
## Chain 2 Iteration:  400 / 2000 [ 20%]  (Warmup) 
## Chain 2 Iteration:  500 / 2000 [ 25%]  (Warmup) 
## Chain 2 Iteration:  600 / 2000 [ 30%]  (Warmup) 
## Chain 2 Iteration:  700 / 2000 [ 35%]  (Warmup) 
## Chain 2 Iteration:  800 / 2000 [ 40%]  (Warmup) 
## Chain 2 Iteration:  900 / 2000 [ 45%]  (Warmup) 
## Chain 2 Iteration: 1000 / 2000 [ 50%]  (Warmup) 
## Chain 2 Iteration: 1001 / 2000 [ 50%]  (Sampling) 
## Chain 2 Iteration: 1100 / 2000 [ 55%]  (Sampling) 
## Chain 2 Iteration: 1200 / 2000 [ 60%]  (Sampling) 
## Chain 2 Iteration: 1300 / 2000 [ 65%]  (Sampling) 
## Chain 2 Iteration: 1400 / 2000 [ 70%]  (Sampling) 
## Chain 2 Iteration: 1500 / 2000 [ 75%]  (Sampling) 
## Chain 2 Iteration: 1600 / 2000 [ 80%]  (Sampling) 
## Chain 2 Iteration: 1700 / 2000 [ 85%]  (Sampling) 
## Chain 2 Iteration: 1800 / 2000 [ 90%]  (Sampling) 
## Chain 2 Iteration: 1900 / 2000 [ 95%]  (Sampling) 
## Chain 2 Iteration: 2000 / 2000 [100%]  (Sampling) 
## Chain 2 finished in 0.1 seconds.
## Chain 3 Iteration:    1 / 2000 [  0%]  (Warmup) 
## Chain 3 Iteration:  100 / 2000 [  5%]  (Warmup) 
## Chain 3 Iteration:  200 / 2000 [ 10%]  (Warmup) 
## Chain 3 Iteration:  300 / 2000 [ 15%]  (Warmup) 
## Chain 3 Iteration:  400 / 2000 [ 20%]  (Warmup) 
## Chain 3 Iteration:  500 / 2000 [ 25%]  (Warmup) 
## Chain 3 Iteration:  600 / 2000 [ 30%]  (Warmup) 
## Chain 3 Iteration:  700 / 2000 [ 35%]  (Warmup) 
## Chain 3 Iteration:  800 / 2000 [ 40%]  (Warmup) 
## Chain 3 Iteration:  900 / 2000 [ 45%]  (Warmup) 
## Chain 3 Iteration: 1000 / 2000 [ 50%]  (Warmup) 
## Chain 3 Iteration: 1001 / 2000 [ 50%]  (Sampling) 
## Chain 3 Iteration: 1100 / 2000 [ 55%]  (Sampling) 
## Chain 3 Iteration: 1200 / 2000 [ 60%]  (Sampling) 
## Chain 3 Iteration: 1300 / 2000 [ 65%]  (Sampling) 
## Chain 3 Iteration: 1400 / 2000 [ 70%]  (Sampling) 
## Chain 3 Iteration: 1500 / 2000 [ 75%]  (Sampling) 
## Chain 3 Iteration: 1600 / 2000 [ 80%]  (Sampling) 
## Chain 3 Iteration: 1700 / 2000 [ 85%]  (Sampling) 
## Chain 3 Iteration: 1800 / 2000 [ 90%]  (Sampling) 
## Chain 3 Iteration: 1900 / 2000 [ 95%]  (Sampling) 
## Chain 3 Iteration: 2000 / 2000 [100%]  (Sampling) 
## Chain 3 finished in 0.1 seconds.
## Chain 4 Iteration:    1 / 2000 [  0%]  (Warmup) 
## Chain 4 Iteration:  100 / 2000 [  5%]  (Warmup) 
## Chain 4 Iteration:  200 / 2000 [ 10%]  (Warmup) 
## Chain 4 Iteration:  300 / 2000 [ 15%]  (Warmup) 
## Chain 4 Iteration:  400 / 2000 [ 20%]  (Warmup) 
## Chain 4 Iteration:  500 / 2000 [ 25%]  (Warmup) 
## Chain 4 Iteration:  600 / 2000 [ 30%]  (Warmup) 
## Chain 4 Iteration:  700 / 2000 [ 35%]  (Warmup) 
## Chain 4 Iteration:  800 / 2000 [ 40%]  (Warmup) 
## Chain 4 Iteration:  900 / 2000 [ 45%]  (Warmup) 
## Chain 4 Iteration: 1000 / 2000 [ 50%]  (Warmup) 
## Chain 4 Iteration: 1001 / 2000 [ 50%]  (Sampling) 
## Chain 4 Iteration: 1100 / 2000 [ 55%]  (Sampling) 
## Chain 4 Iteration: 1200 / 2000 [ 60%]  (Sampling) 
## Chain 4 Iteration: 1300 / 2000 [ 65%]  (Sampling) 
## Chain 4 Iteration: 1400 / 2000 [ 70%]  (Sampling) 
## Chain 4 Iteration: 1500 / 2000 [ 75%]  (Sampling) 
## Chain 4 Iteration: 1600 / 2000 [ 80%]  (Sampling) 
## Chain 4 Iteration: 1700 / 2000 [ 85%]  (Sampling) 
## Chain 4 Iteration: 1800 / 2000 [ 90%]  (Sampling) 
## Chain 4 Iteration: 1900 / 2000 [ 95%]  (Sampling) 
## Chain 4 Iteration: 2000 / 2000 [100%]  (Sampling) 
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.9 seconds.
 summary(dataset1.brm.bayesiant)
##  Family: student 
##   Links: mu = identity; sigma = identity; nu = identity 
## Formula: effectiveness ~ condition 
##    Data: exp1.data (Number of observations: 123) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Population-Level Effects: 
##                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept             6.63      0.18     6.29     6.97 1.00     3579     2545
## conditionno_graph     0.16      0.24    -0.31     0.63 1.00     3592     2585
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     1.29      0.12     1.07     1.52 1.00     3264     2435
## nu       16.53     11.57     4.27    46.43 1.00     3170     2975
## 
## Samples were drawn using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
color_scheme_set(scheme = "purple")
plot(dataset1.brm.bayesiant)

  dataset1.bayesiant.y <- exp1.data$effectiveness
  dataset1.bayesiant.yrep <- posterior_predict(dataset1.brm.bayesiant, draws = 50)
       
  kable(head(as_tibble(dataset1.bayesiant.yrep[,1:11])))            
## Warning: The `x` argument of `as_tibble.matrix()` must have unique column names if `.name_repair` is omitted as of tibble 2.0.0.
## Using compatibility `.name_repair`.
V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11
7.350071 7.550894 4.346016 7.402838 6.105783 6.529116 8.815189 8.225576 6.937635 7.539691 5.342559
7.294639 7.101592 4.382338 8.122888 9.732572 5.105547 5.986431 7.017659 7.960890 6.007279 5.592948
6.181452 2.912017 5.652890 6.575436 7.737640 7.520308 8.458674 6.613323 8.191181 7.433597 6.322437
8.535437 7.939905 5.246086 5.892014 4.370619 6.351047 6.903111 4.641686 7.422193 7.811355 5.179726
7.361354 4.894749 6.116840 5.985281 6.603010 6.110403 6.298743 7.107104 5.918865 6.950818 6.126804
6.634058 8.767211 7.708752 4.349791 6.250896 5.995275 6.298566 9.387987 5.926242 6.926282 8.637597
ppc_hist(y = dataset1.bayesiant.y,
         yrep = dataset1.bayesiant.yrep[1:8, ],
         binwidth = .5)

ppc_dens_overlay(y = dataset1.bayesiant.y,
                 yrep = dataset1.bayesiant.yrep[1:30, ])

ppc_stat_grouped(y = dataset1.bayesiant.y,
                 yrep = dataset1.bayesiant.yrep,
                 group = exp1.data$condition)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Compared to the frequentist estimates:

bind_rows(tidy.wrapper(dataset1.lm.freqt),
          tidy.wrapper(dataset1.brm.bayesiant)) %>%
  ggplot() +
  geom_pointrange(
    aes(
      x = model,
      y = estimate,
      ymin = conf.low,
      ymax = conf.high,
      color = term
    ),
    position = position_dodge(width = 0.2)
  ) +
  scale_color_brewer(palette = "Set1") +
  ylab('effectiveness') +
  scale_y_continuous(breaks = 1:9, limits = c(1, 9)) +
  coord_flip()

dataset1.bayesiant.posterior_epred <- tibble(condition = c('graph', 'no_graph')) %>%
  add_epred_draws(dataset1.brm.bayesiant, re_formula = NA, allow_new_levels = FALSE) %>%
  ungroup()

kable(head(dataset1.bayesiant.posterior_epred))
condition .row .chain .iteration .draw .epred
graph 1 NA NA 1 6.50377
graph 1 NA NA 2 6.64234
graph 1 NA NA 3 6.59090
graph 1 NA NA 4 6.50090
graph 1 NA NA 5 6.69152
graph 1 NA NA 6 6.47357
dataset1.bayesiant.posterior_comparison <- dataset1.bayesiant.posterior_epred %>%
  select(-c(.chain, .iteration, .row)) %>% 
  compare_levels(variable = .epred, by = condition)

kable(head(dataset1.bayesiant.posterior_comparison))
.draw condition .epred
1 no_graph - graph 0.1420210
2 no_graph - graph 0.1585310
3 no_graph - graph 0.3108710
4 no_graph - graph 0.0574165
5 no_graph - graph 0.3252280
6 no_graph - graph 0.2323100
kable(dataset1.bayesiant.posterior_comparison %>%
           mean_qi(.epred))
condition .epred .lower .upper .width .point .interval
no_graph - graph 0.1609358 -0.3058675 0.6307221 0.95 mean qi
dataset1.bayesiant.posterior_comparison %>%
  median_qi(.epred) %>%
  ggplot() +
  geom_point(aes(x = .epred, y = condition), color = DARK_PURPLE, size = 3) +
  geom_errorbarh(
    aes(xmin = .lower, xmax = .upper, y = condition),
    color = DARK_PURPLE,
    alpha = .5,
    size = 2,
    height = 0
  ) +
  geom_vline(aes(xintercept = 0), linetype = "dashed", color = "gray") +
  coord_cartesian(ylim = c(0, 2), xlim = c(-1, 1))  +
  xlab('') + ylab('')

Model 3. Ordinal Logistic Regression

dataset1.brm.olr.formula <- 
   bf(
    effectiveness ~ condition,
    family = cumulative("logit")
  )
get_prior(dataset1.brm.olr.formula, data = exp1.data)
##                 prior     class              coef group resp dpar nlpar bound
##                (flat)         b                                              
##                (flat)         b conditionno_graph                            
##  student_t(3, 0, 2.5) Intercept                                              
##  student_t(3, 0, 2.5) Intercept                 1                            
##  student_t(3, 0, 2.5) Intercept                 2                            
##  student_t(3, 0, 2.5) Intercept                 3                            
##  student_t(3, 0, 2.5) Intercept                 4                            
##  student_t(3, 0, 2.5) Intercept                 5                            
##  student_t(3, 0, 2.5) Intercept                 6                            
##  student_t(3, 0, 2.5) Intercept                 7                            
##  student_t(3, 0, 2.5) Intercept                 8                            
##        source
##       default
##  (vectorized)
##       default
##  (vectorized)
##  (vectorized)
##  (vectorized)
##  (vectorized)
##  (vectorized)
##  (vectorized)
##  (vectorized)
##  (vectorized)
dataset1.brm.olr.priors = c(
      prior(normal(0, 1), class = "b"),
      prior(student_t(3, 0, 1), class = "Intercept")
)
dataset1.brm.olr.priorchecks <- brm(
    effectiveness ~ condition,
    prior = dataset1.brm.olr.priors,
    data = exp1.data,
    family= cumulative("logit"),
    backend = 'cmdstanr',
    sample_prior = 'only',
    file = "02-bayesian_stats/rds/dataset1.brm.bayesiant.priorchecks.rds"
    )



dataset1.olr.yprior <-
  posterior_predict(dataset1.brm.olr.priorchecks)


ggplot() +
  geom_histogram(aes(x = dataset1.olr.yprior),
               fill = '#351c75',
               alpha = .5,
               size = 1) +
  scale_x_continuous(breaks = 1:9, limits = c(1,9)) + 
  xlab('prior draws') +
  ggtitle('prior predictive check')
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 165719 rows containing non-finite values (stat_bin).
## Warning: Removed 2 rows containing missing values (geom_bar).

dataset1.brm.olr1 = 
  brm(dataset1.brm.olr.formula,
    prior = dataset1.brm.olr.priors,
    data = exp1.data,
    backend = "cmdstanr",
    file = "02-bayesian_stats/rds/dataset1.brm.olr1.rds"
  )

summary(dataset1.brm.olr1)
## Warning: There were 27 divergent transitions after warmup. Increasing
## adapt_delta above may help. See http://mc-stan.org/misc/warnings.html#divergent-
## transitions-after-warmup
##  Family: cumulative 
##   Links: mu = logit; disc = identity 
## Formula: effectiveness ~ condition 
##    Data: exp1.data (Number of observations: 123) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Population-Level Effects: 
##                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept[1]         -4.56      0.96    -6.78    -2.99 1.00     1070      513
## Intercept[2]         -3.99      0.78    -5.82    -2.71 1.00     1583     1198
## Intercept[3]         -3.29      0.56    -4.55    -2.32 1.00     2496     2468
## Intercept[4]         -2.14      0.34    -2.85    -1.50 1.00     3348     3016
## Intercept[5]         -1.26      0.26    -1.79    -0.76 1.00     3758     3330
## Intercept[6]         -0.27      0.23    -0.74     0.18 1.00     3652     3366
## Intercept[7]          1.26      0.26     0.77     1.81 1.00     3532     3334
## Intercept[8]          2.22      0.33     1.62     2.91 1.00     3937     2992
## conditionno_graph     0.18      0.31    -0.40     0.81 1.00     3434     2673
## 
## Family Specific Parameters: 
##      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc     1.00      0.00     1.00     1.00 1.00     4000     4000
## 
## Samples were drawn using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
dataset1.brm.olr2 = 
  brm(dataset1.brm.olr.formula,
    prior = dataset1.brm.olr.priors,
    data = exp1.data,
    backend = "cmdstanr",
    warmup = 1500,
    iter = 2500,
    control = list(adapt_delta = 0.99, max_treedepth = 15),
    file = "02-bayesian_stats/rds/dataset1.brm.olr2.rds"
  )

summary(dataset1.brm.olr2)
##  Family: cumulative 
##   Links: mu = logit; disc = identity 
## Formula: effectiveness ~ condition 
##    Data: exp1.data (Number of observations: 123) 
## Samples: 4 chains, each with iter = 2500; warmup = 1500; thin = 1;
##          total post-warmup samples = 4000
## 
## Population-Level Effects: 
##                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept[1]         -4.58      1.00    -6.88    -3.00 1.00     2194     1948
## Intercept[2]         -3.99      0.78    -5.76    -2.71 1.00     2806     2384
## Intercept[3]         -3.30      0.55    -4.48    -2.32 1.00     3705     2806
## Intercept[4]         -2.15      0.34    -2.85    -1.51 1.00     4378     3270
## Intercept[5]         -1.26      0.26    -1.78    -0.76 1.00     4525     3093
## Intercept[6]         -0.27      0.24    -0.73     0.19 1.00     3935     3132
## Intercept[7]          1.27      0.26     0.76     1.79 1.00     4089     3216
## Intercept[8]          2.23      0.34     1.58     2.92 1.00     4890     3122
## conditionno_graph     0.19      0.31    -0.43     0.82 1.00     4856     2853
## 
## Family Specific Parameters: 
##      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc     1.00      0.00     1.00     1.00 1.00     4000     4000
## 
## Samples were drawn using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
plot(dataset1.brm.olr2)

dataset1.olr.y <- exp1.data$effectiveness
dataset1.olr.yrep <- posterior_predict(dataset1.brm.olr2)

ppc_hist(y = dataset1.olr.y,
         yrep = dataset1.olr.yrep[1000:1007,], binwidth = .5)

ppc_dens_overlay(y = dataset1.olr.y,
         yrep = dataset1.olr.yrep[2000:2030,])

ppc_stat_grouped(y = dataset1.olr.y,
         yrep = dataset1.olr.yrep, group = exp1.data$condition)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

dataset1.olr.posterior_epred <-
  tibble(condition = c('graph', 'no_graph')) %>%
  add_epred_draws(dataset1.brm.olr2,
                  re_formula = NA,
                  allow_new_levels = FALSE) %>%
  ungroup()

kable(head(dataset1.olr.posterior_epred))
condition .row .chain .iteration .draw .category .epred
graph 1 NA NA 1 1 0.0045580
graph 1 NA NA 2 1 0.0063148
graph 1 NA NA 3 1 0.0069149
graph 1 NA NA 4 1 0.0391315
graph 1 NA NA 5 1 0.0081217
graph 1 NA NA 6 1 0.0248435
dataset1.olr.posterior_comparison <-
  dataset1.olr.posterior_epred %>%
  select(-c(.chain, .iteration, .row)) %>%
  group_by(.category) %>%
  compare_levels(variable = .epred, by = condition)

kable(head(dataset1.olr.posterior_comparison %>%
             mean_qi()))
.category condition .epred .lower .upper .width .point .interval
1 no_graph - graph -0.0025112 -0.0174983 0.0069861 0.95 mean qi
2 no_graph - graph -0.0013803 -0.0110381 0.0040049 0.95 mean qi
3 no_graph - graph -0.0026421 -0.0162269 0.0074813 0.95 mean qi
4 no_graph - graph -0.0102726 -0.0506452 0.0227708 0.95 mean qi
5 no_graph - graph -0.0137966 -0.0621373 0.0329472 0.95 mean qi
6 no_graph - graph -0.0145299 -0.0680404 0.0330627 0.95 mean qi
dataset1.olr.posterior_comparison %>%
  mean_qi(.epred) %>%
  ggplot() +
  geom_point(aes(y = .epred, x = .category), size = 3, color = DARK_PURPLE) +
  geom_errorbar(
    aes(ymin = .lower, ymax = .upper, x = .category),
    width = 0,
    size = 2,
    color = DARK_PURPLE,
    alpha = .5
  ) +
  geom_hline(aes(yintercept = 0), linetype = "dashed", color = "gray") +
  xlab('') + ylab('no_graph - graph') +
  theme(axis.title.y = element_text(angle = 0, vjust = .5))

# Dataset 2

Load the data

dataset2 = readr::read_csv("02-bayesian_stats/data/exp2.csv") %>%
  mutate(condition = condition == 'expansive') %>%
  group_by(participant)
## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   participant = col_double(),
##   adjP10 = col_double(),
##   adjP20 = col_double(),
##   adjP30 = col_double(),
##   adjPDiff = col_double(),
##   condition = col_character(),
##   gender = col_character(),
##   index = col_double(),
##   change = col_double(),
##   subgroup = col_character(),
##   orig = col_double(),
##   BIS = col_character()
## )
kable(head(dataset2))
participant adjP10 adjP20 adjP30 adjPDiff condition gender index change subgroup orig BIS
1 56.50000 55.80000 58.40000 1.90000 TRUE female 69 3.362832 expansiveHigh 56.90000 high
2 27.87500 40.14286 36.00000 8.12500 FALSE female 60 29.147982 constrictiveLow 34.67262 low
3 61.00000 59.00000 76.50000 15.50000 TRUE female 55 25.409836 expansiveLow 65.50000 low
4 24.57143 32.75000 37.85714 13.28571 FALSE female 79 54.069767 constrictiveHigh 31.72619 high
5 64.71429 54.00000 41.00000 -23.71429 TRUE female 69 -36.644592 expansiveHigh 53.23810 high
6 35.14286 52.40000 45.60000 10.45714 FALSE female 43 29.756098 constrictiveLow 44.38095 low
dataset2 %>%
  mutate(c = as.factor(condition)) %>%
  ggplot(aes(x = change)) +
  geom_histogram(aes(y = ..density..),
                 binwidth = 6,
                 fill = DARK_PURPLE,
                 color = 'white') +
  geom_density(size = 1, adjust = 1.5, color = 'lightgray') +
  scale_x_continuous(limits = c(-200, 200)) 
## Warning: Removed 1 rows containing non-finite values (stat_bin).
## Warning: Removed 1 rows containing non-finite values (stat_density).
## Warning: Removed 2 rows containing missing values (geom_bar).

Model 1

The first model is the BEST test model as described by Kruschke in the paper Bayesian estimation supersedes the t-test. In this model, \(\beta\) indicates the mean difference in the outcome variable between the two groups (in this case, the percent change in the BART scores). We fit different priors on \(\beta\) and set different weights on these priors to obtain our posterior estimate.

\[ \begin{align} y_{i} &\sim \mathrm{T}(\nu, \mu, \sigma) \\ \mu &= \alpha_{0} + \beta * x_i \\ \sigma &= \sigma_{a} + \sigma_{b}*x_i \\ \beta &\sim \mathrm{N}(\mu_{0}, \sigma_{0}) \\ \sigma_a, \sigma_b &\sim \mathrm{Cauchy}(0, 2) \\ \nu &\sim \mathrm{exp}(1/30) \end{align} \]

dataset2.brm.student.formula <- bf(
    change ~ condition, 
    sigma ~ condition, 
    family = student()
  )

kable(head(as_tibble(get_prior(dataset2.brm.student.formula, data = dataset2))))
prior class coef group resp dpar nlpar bound source
b default
b conditionTRUE default
student_t(3, 22.6, 38.8) Intercept default
gamma(2, 0.1) nu default
b sigma default
b conditionTRUE sigma default
dataset2.brm.student.priorchecks = brm(
  dataset2.brm.student.formula,
  prior = c(
    prior(normal(0, 2), class = "b"),
    prior(cauchy(0, 2), class = "b", dpar = "sigma"),
    prior(exponential(0.033), class = "nu"),
    prior(student_t(3, 0, 5), class = "Intercept"),
    prior(student_t(3, 0, 2), class = "Intercept", dpar = "sigma")
  ),
  data = dataset2,
  backend = "cmdstanr",
  sample_prior = 'only',
  file = "02-bayesian_stats/rds/dataset2.brm.student.priorchecks.rds"
)

# dataset2.student.yprior <-
#   posterior_predict(dataset2.brm.student.priorchecks)

# ggplot() +
#   geom_density(
#     aes(x = dataset2.student.yprior),
#     color = '#351c75',
#     alpha = .5,
#     size = 1
#   ) +
#   xlab('prior draws') +
#   ggtitle('prior preditive check')
dataset2.brm.student = brm(
 bf(change ~ condition, 
 sigma ~ condition,
 family = student()),
 prior = c(
   prior(normal(0, 2), class = "b"),
   prior(cauchy(0, 2), class = "b", dpar = "sigma"),
   prior(exponential(.033), class = "nu"),
   prior(student_t(3, 0, 5), class = "Intercept"),
   prior(student_t(3, 0, 2), class = "Intercept", dpar = "sigma")
 ), 
 data = dataset2,
 backend = "cmdstanr",
 file = "02-bayesian_stats/rds/dataset2.brm.student.rds"
)

dataset2.brm.student$prior
##               prior     class          coef group resp  dpar nlpar bound
##        normal(0, 2)         b                                           
##        normal(0, 2)         b conditionTRUE                             
##        cauchy(0, 2)         b                          sigma            
##        cauchy(0, 2)         b conditionTRUE            sigma            
##  student_t(3, 0, 5) Intercept                                           
##  student_t(3, 0, 2) Intercept                          sigma            
##  exponential(0.033)        nu                                           
##        source
##          user
##  (vectorized)
##          user
##  (vectorized)
##          user
##          user
##          user
summary(dataset2.brm.student)
##  Family: student 
##   Links: mu = identity; sigma = log; nu = identity 
## Formula: change ~ condition 
##          sigma ~ condition
##    Data: dataset2 (Number of observations: 80) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Population-Level Effects: 
##                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept              21.15      4.95    11.41    31.04 1.00     2751     2071
## sigma_Intercept         3.37      0.20     2.98     3.76 1.00     2331     2293
## conditionTRUE          -0.15      1.92    -4.01     3.46 1.00     3711     2961
## sigma_conditionTRUE     0.14      0.22    -0.31     0.55 1.00     3648     2440
## 
## Family Specific Parameters: 
##    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## nu     4.87      5.87     1.59    16.82 1.00     1781     1303
## 
## Samples were drawn using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
plot(dataset2.brm.student)

dataset2.student.y <- dataset2.brm.student$data$change
dataset2.student.yrep <- posterior_predict(dataset2.brm.student)

dataset2.student.epred <- tibble(condition = c(TRUE, FALSE)) %>%
  add_epred_draws(dataset2.brm.student, re_formula = NA, allow_new_levels = TRUE) %>%
  ungroup()
ppc_hist(y = dataset2.student.y,
         yrep = dataset2.student.yrep[100:107,], binwidth = 10)

ppc_dens_overlay(y = dataset2.student.y,
         yrep = dataset2.student.yrep[2000:2030,])

ppc_stat_grouped(y = dataset2.student.y,
         yrep = dataset2.student.yrep,
         group = dataset2$condition,
         binwidth = 5)

dataset2.student.epred_comparison <- tibble(condition = c(TRUE, FALSE)) %>%
 add_epred_draws(dataset2.brm.student, re_formula = NA, allow_new_levels = FALSE) %>%
 ungroup() %>%
 select(-c(.chain, .iteration, .row)) %>% 
 compare_levels(variable = .epred, by = condition) %>%
 rename(diff = .epred)


kable(head(dataset2.student.epred_comparison))
.draw condition diff
1 TRUE - FALSE -2.918310
2 TRUE - FALSE 2.367050
3 TRUE - FALSE 3.045840
4 TRUE - FALSE -0.455904
5 TRUE - FALSE -1.398820
6 TRUE - FALSE 2.193080
dataset2.student.epred_comparison %>%
  select(diff) %>%
  mean_qi() %>%
ggplot() +
  geom_point(aes(x = diff, y = condition), size = 3, color = DARK_PURPLE) + 
  geom_errorbarh(
    aes(xmin = .lower, xmax = .upper, y = condition),
    height = 0,
    color = DARK_PURPLE,
    alpha = .5,
    size = 2
  ) +
  geom_vline(aes(xintercept = 0), linetype = "dashed", color = "gray") + 
  coord_cartesian(ylim = c(0, 2), xlim = c(-5, 5))  + 
  scale_y_discrete(label = c('expansive - not expansive')) + 
  xlab('')  + ylab('')
## Adding missing grouping variables: `condition`